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Abstract 

The classical problem of two uniformly charged infinite planes in electrochemical equilibrium with 
an infinite monovalent salt reservoir is solved exactly at the mean-field nonlinear Poisson-Boltzmann 
(PB) level, including an explicit expression of the associated nonlinear electrostatic contribution to the 
semi-grand-canonical potential. A linearization of the nonlinear functional is presented that leads to 
Debye-Hiickel-like equations agreeing asymptotically with the nonlinear PB results in the weak-coupling 
(high-temperature) and counterionic ideal-gas limits. This linearization scheme yields artifacts in the low- 
temperature, large-separation or high-surface charge limits. In particular, the osmotic-pressure difference 
between the interplane region and the salt reservoir becomes negative in the above limits, in disagreement 
with the exact (at mean-field level) nonlinear PB solution. By using explicitly gauge-invariant forms of 
the electrostatic potential we show that these artifacts — although thermodynamically consistent with 
quadratic expansions of the nonlinear functional — can be traced back to the non-fulfillment of the 
underlying assumptions of the linearization. Explicit comparison between the analytical expressions 
of the exact nonlinear solution and the corresponding linearized equations allows us to show that the 
linearized results are asymptotically exact in the weak-coupling and counterionic ideal-gas limits, but 
always fail otherwise, predicting negative osmotic-pressure differences. By taking appropriate limits of the 
full nonlinear PB solution, we provide asymptotic expressions for the semi-grand-canonical potential and 
the osmotic-pressure difference that involve only elementary functions, which cover the complementary 
region where the linearized theory breaks down. 

1 Introduction 

In the preceding paperS we performed a linearization of the mean-field Poisson-Boltzmann (PB) density 
functional for spherical Wigner-Seitz cells that agrees asymptotically with the PB results in the weak- 
coupling (high-temperature) limit by adopting an explicitly gauge-invariant approach. However, as already 
pointed out previously in the literature,Btj the linearization scheme yields artifacts in the low-temnerajbjre, 
high-surface charge or infinite-dilution (of polyions) limits for the Donnan equilibrium problemoQ&Lrlj in 
spherical geometry, which describes a suspension of spherical charged polyions in electrochemical equilibrium 
with an infinite salt reservoir. In these limits the linearized osmotic-pressure difference between the colloidal 
suspension and the salt reservoir becomes negative, in disagreement with the full nonlinear PB result, that 
always displays positive osmotic-pressure differences. Because the nonlinear PB equation is not analytically 
solvable in spherical geometry — even in the simplest salt-free case, when only neutralizing counterions are 
present — one must rely on numerical calculations to establish comparisons between the nonlinear and the 
linearized equations. This motivated us to consider the prototype case represented by the classical problem 
of two uniformly charged infinite planes in electrochemical equilibrium with an infinite salt reservoir, when 
the exact analytical solution of the nonlinear problem is possible. The explicit analytical comparison between 
the exact (at the mean-field level) full nonlinear and the approximated linearized equations allows us to trace 
back the underlying reasons of the breakdown of the linearization scheme that is intrinsically associated with 
its range of validity. Additionally, the study of this exactly-solvable case sheds some light on the question of 
the proper definition of the linearized osmotic pressure that was previously considered in Ref. |3| . 

Moreover, to our knowledge, the explicit calculation of the semi-grand-canonical potential for two uni- 
formly charged infinite planes at the nonlinear mean-field PB level has only been reported in connection 
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to the polyelectrolyte-brush problem.EI In that work, however, the thermodynamical potential also included 
electrostatic and elastic contributions arising from the polyelectrolyte brushes, and therefore, these need to 
be subtracted out. The knowledge of the thermodynamic potential allows us to derive all thermodynamic 
properties of the two charged infinite plpins-jppblcm. We note that it can be also extended to curved surfaces 
by using the Derjaguin approximationil3'tiH13 It is then possible to determine the normal forces (per unit 
area) between these surfaces when their separation distance is much smaller than their curvature radius. We 
will present the exact nonlinear semi-grand-canonical functional from which we derive approximate expres- 
sions. These involve only elementary functions and provide excellent approximations to the full nonlinear 
PB results within the whole range of parameters. 

The remainder of the paper is organized as follows. In Section |^ the model is introduced and the 
exact nonlinear solution is obtained. In Section || the linearization of the appropriate semi-grand-canonical 
functional is performed. In Section |4] we present explicit analytic comparisons between the exact nonlinear 
and the linearized equations of state. In Section |^ we give some concluding remarks. In Appendix ^| 
we compare the exact nonlinear and the self-consistent linearized averaged densities used to perform the 
quadratic expansion of the nonlinear functional. In Appendix [b| we derive the asymptotic expansions of 
the exact nonlinear solution. In Appendix |c| we obtain extended expansions of the exact nonlinear solution 
that are valid in the Gouy-Chapman (high-surface charge) and the large-separation limits, regions where the 
linearization scheme breaks down. In Appendix [S] we present linearized equations that preserve the exactness 
of the Legendre transformation between the semi-grand-canonical and canonical-ensemble formulations of 
the problem. 



2 Exact nonlinear Poisson-Boltzmann solution 

The system to be considered is comprised of two infinite planar surfaces a distance 2L apart, each with 
surface charge — aq, where q > is the elementary charge, in electrochemical equilibrium with an infinite 
monovalent salt reservoir of bulk salt density n^. The microions (positive counterions and salt ions) are free 
to move in the region —L < x < L between the two charged surfaces, where we introduced a Cartesian axis 
perpendicular to the planes in which the midplanc is located at x — and the two charged planes at x = ±L. 
At the mean-field level of approximation the ions are treated as inhomogeneous ideal gases described by their 
average local number densities n±(x). We do not distinguish between (positive) counterions and positive 
ions derived from the salt dissociation. The total charge number density (counterions, salt ions and the 
negative surface charge on the planes) of the system, 

p(x) = n + (x) — ri— (x) — aS(x + L) — a8{x — L), (1) 

where 6 is the one-dimensional Dirac delta- function, is related to the reduced electrostatic potential ip(x) = 
(3q^(x), which satisfies the (exact) Poisson equation, 

^=-^BP(x), (2) 

where £b = Pq 2 A i s the Bjcrrum length and = fc^T is the thermal energy. It is implicitly assumed that 
the solvent dielectric constant e remains the same outside the region containing the salt solution (|x| > L), 
so image-charge effects due to dielectric contrast are absent. The formal solution of the Poisson equation (||) 
may be written in terms of the one-dimensional Green function G\(x, xo). 
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ip(x) = r lim Ib I dx Gi(x,x ) p(x ), — j = -4ir5(x - x ), (3) 



which in turn allows us to express the mean-field semi-grand-canonical functional (for one charged plane) 
per unit area as 
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where C± are the thermal de Broglie wavelengths of cations (including the positive counterions) and anions, 
respectively, and the (mean-field) microion chemical potentials /3[i± = ln(ribC±) assume ideal gases of 
uniform density n-^ for both types of ions in the infinite salt reservoir. The last equality results from 
integrating the electrostatic energy (per unit area) — the first term of Eq. (||) — by parts and using the fact 
that surface contributions vanish due to Gauss' law and the overall electroneutrality of the system, 
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lim / da; p(x) = 0, 



or 



da; [n+(x) — ri-(x)] = a. 



(5) 



The equilibrium density profiles are obtained by minimizing the PB semi-grand-canonical functional (|]) 
under the charge- neutrality constraint (0) , 
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where we introduced the Lagrange multiplier p c \ to ensure the neutrality condition (^|). This yields the 
Boltzmann- weighted ionic profiles, 



n± (x) = n h exp [±/?/i cl T ip(x)] , 
and the Lagrange multiplier fj, e \ is found by imposing the charge- neutrality condition (|^) , 



yjnl + (2?i b ) 2 a+a- ± n c ± 



where n c = a/L is the average density of counterions in the interplane region 
denote unweighted spatial averages over the volume available to the microions, 
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< L and the brackets 
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In close analogy to the spherical case — cf. Appendix E of Ref. jy — it is possible to write the nonlinear 
equilibrium density profiles in an explicitly gauge-invariant form by inserting the Lagrange multiplier (ph 
into the Boltzmann- weighted ionic profiles, Eqs. (0), 
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By explicitly gauge-invariant we mean that the equilibrium profiles do not depend explicitly on a particular 
choice of the zero of the potential, because they depend only on the difference (-0) — Henceforth 
gauge-invariant will be a short writing to explicitly gauge-invariant. In particular, in the salt-free (rib — > 0) 
limit, these gauge-invariant forms lead — in a direct and transparent way — to the vanishing coion profile, 
ri-(x) = 0, and to the salt-free equilihriiim counterion profile, n + (x) — n c e^~^^ /a + . 

The most commonly used gaugeE3'E3 is the one in which the charge-neutrality Lagrange multiplier is 
zero, /i c i = 0, which does not correspond to the gauge in which the electrostatic potential at the infinite 
salt reservoir vanishesO Henceforth, to simplify the notation and the calculations we will use the standard 
gauge fi e \ = to treat the nonlinear problem. We should keep in mind, however, that the fixed-gauge 
electrostatic potential ip(x) = ijj{x) — (3p c \ will no longer be gauge- invariant: its value at a particular point 
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- let us say, at the midplane, tpo = ip(x — 0), or at the charged surfaces, pi, = ip(x — L) — will be 
determined by imposing the overall charge neutrality (|^) in the whole system. They can no longer be chosen 
arbitrarily, in contrast to their gauge-invariant counterparts, f/'o = ^(x = 0) or ipL = = L). In the 
gauge-invariant formulation, either ipo or ipi, may be chosen arbitrarily — but not both simultaneously — 
because the difference ipL ~ V'o ~ P>l ~ <A) must eventually be preserved. 

In the standard, gauge the nonlinear problem reduces into solving the usual PB equation for two charged 
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with the appropriate boundary conditions 



k 2 sinh ip{x) + — [S(x + L) + S(x - L)] , 
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ip'(x = 0) = ip'(x = 0) = 0, 



and 



if/(x = TL) = <p'(x = TL) = ±j, 



(13) 



the prime (') denoting differentiation with respect to the argument. We have introduced two length scales: 
the Debye screening length associated with the bulk density n b of the infinite salt reservoir, 
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which gives the characteristic (algebraic) decay length of the counterion distribution (for a salt-free system 
around an infinite charged plane with bare surface charge a. 
Using the mathematical identity 
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dx 2 2 dip 

it is possible to integrate the nonlinear PB Eq. (|lj) exactly, 
[<//(a;)] 2 = k 2 [2cosh<^(ir) — 2 cosh po] , 
dip 
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whose solution is written in terms of the midplane electrostatic potential < and F(ip\m) 



Jq d6 / \J 1 — m sin 2 9 is the incomplete elliptic integral of the first kind-llS'llS'SSEil Applying the boundary 
conditions (|l^) yields 
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where we defined the two dimensionless distances A = K b A and I = K b L, and p>L < Po < is the surface 
electrostatic potential at the charged planes. Introducing the variable 



i^sinh 2 ^, 



(21) 



which a posteriori will be identified with half of the (dimensionless) osmotic-pressure difference between the 
interplane region and the salt reservoir — cf. Eq. (|33| ) — the two boundary conditions can be combined into 
the eigenvalue equation, 
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where cs(ii 
functions Jl* 



= r £n(u\m) / sn(u|m) is the ratio of the cosine-amplitude and sine-amplitude Jacobi elliptic 
E2rEd The explicit exact solution of the nonlinear PB problem can then be written as 



f(x) — — 2 arcsinh 



Vt 



K h \x\^/1 + t 



1 + t 



\x\ < L. 



(23) 



It should be remarked that th&jexact solution to the-nonlinear PB problem may be cast in several equivalent 
forms. Verwey and Overbeek,E3 quoted by Hunter ,E3 gave an alternative form for the implicit solution 
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where K(m) = FLjl/2 \ to) is the complete elliptic integral of the first kind' 16 
Borkovec's versiorO to the explicit solution (E31) reads 



-2tp 

(24) 

11 while Behrens and 



cp(x) = (^ + 21ncd ( C - ipo/2 K h \x\/2 \ e 2vo ) 



(25) 

where cd(u|m) is the cd Jacobi elliptic function.H3ll3E3'Ell However, none of these previous works presented 
the explicit expression for the nonlinear PB semi-grand-canonical potential f2 = n[n±(x)] cqu n, which can 
be extracted from Ref. Q by neglecting the electrostatic and elastic contributions arising from the polyelec- 
trolyte brushes. ._. 

The dimensionless excessEIl semi-grand-canonical potential per unit area, 
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may be evaluated inserting the exact nonlinear solution ( p3| ) into the semi-grand-canonical functional (0) 
and performing the integrations. After some tedious algebra, we obtain 
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where E(ip\m) — d6*\/l — to sin 2 9 is the incomplete elliptic integral of the second kindJl3'H3'E2IEiJ ip^ < 
is the reduced electrostatic surface potential at the charged plane at infinite separation and w(A, oo) represents 
the nonlinear excess self-energy of the system. Using the relations 
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and the fact that <pl < and ifoo < 0, the excess semi-grand-canonical potential uj may be cast in a simpler 
form, 
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(31) 

As already pointed out in the introduction, besides its intrinsic relevance, the semi-grand-canonical potential 
for the planar case has also animportant application: by using the Derjaguin approximationJi3'EilE3 valid 
when the range of interactionsES (of order k^ 1 ) and the separation distance 2L between the two curved 
surfaces are much smaller than their curvature radius a, it is possible to determine their normal force (per 
unit area) T at separation 2L by 
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Returning to the planar case, the osmotic-pressure difference AP between the interplane region and the 
infinite salt reservoir can be written in terms of the midplane reduced electrostatic potential tpo, 
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(33) 



which may be obtained by taking the formal derivative .of ui with respect to I, similarly as performed in 
Appendix A of Ref. [Q. Eq. ( |33| ) is a mean-field versionEZl of the boundary-density theorem, which states 
that the osmotic pressure is simply given by the sum of the microionic densities at the midplane (WS-cell 
boundary). This simple relation does not hold beyond the mean- field level because of finite ionic-size effects 
and the presence of microionic correlations between particles located in the different semi-spaces separated 
by the mid-plane — even though it still does for one charged plane with the electrolyte confined by a neutral 
midplaneO We restrict ourselves, however, to the nonlinear mean- field result (|3q), which clearly predicts 
that the osmotic-pressure difference II is always positive. 

We should remark that the nonlinear osmotic-pressure difference (p3) may also be written in a gauge- 
invariant form using the gauge-invariant equilibrium density profiles (|1 1|) , 
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(34) 



where ipo is the arbitrary midplane electrostatic potential in the gauge-invariant formulation. It does not 
need necessarily to coincide with the midplane electrostatic potential cpo = 2 arcsinh y/IL/2 in the standard 
gauge /i c i = 0. The gauge- invariant form (^) of the nonlinear osmotic-pressure difference will be useful 
later, at the end of Section ^, when establishing a connection between its quadratic expansion about the 
average potential (ip) and its linearized counterpart (54). 



3 Linearization scheme 

To obtain the linearized semi-grand-canonical functional Onn[n±(x)] we truncate the expansion of the PB 
nonlinear semi-grand-canonical functional fl[n±{x)\ 1 Eq. (0), up to second order in the differences n±(x) — 
(n±(x)}, where (n±(x)} = (l/£) J dxn±(x) are the (a priori unknown) average densities, 
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After minimization of the functional I?dh[w±W] with respect to the profiles n±(x) under the overall elec- 
troneutrality constraint (^), [5/5n±(xY\ [J?dh / A — /i c \ J dx p(x)] =0 — analogously as performed for the 
spherical case in Appendix F of Ref. JO — we obtain the self-consistent linearized averaged densities, 



n± = (n±(x)) 1 = 
and the linearized equilibrium density profiles, 



Vn c 2 + (2n b ) 2 ± 



n±(x) =n±[l± (tp(x)) 1 =F ^[x)} , 



(36) 



(37) 



where the notation (• • -) t emphasizes the fact that the self-consistent averaged densities (|36|) were obtained 
within a linearized approximation. Henceforth we will omit the subscript '1' in order to simplify the notation. 
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We should remark that the profile independence of the self-consistent linearized averaged densities given 
by ( |36|) can only be verified after the minimization procedure. When performing the functional minimization 
of the linearized functional J7dh[^±(2;)], one must take the profile dependence of the average expansion 
densities (n± (x)) into account, in addition to the charge-neutrality constraint (||). Although similar quadratic 
expansions about the zero-th order Donnan densities for the planar case were already proposed by Trizac 
and Hansen,E3 they focused their study on finite-size|-effects and did not investigate the consequences of 
the linearization in detail. Deserno and von Grunbergp considered the general d-dimensional problem in a 
fixed-gauge formulation, interpreting these self-consistent linearized averaged densities in terms of an optimal 
linearization point, n± = nbe T ^ opt . 



The linearized expansion densities ( 36 ) , which correspond to the state- independent zero-th order Donnan 



densities, represent the infinite-temperature [Ib = 0) limit of the gauge- invariant forms of the equilibrium 



density profiles (11) and do not coincide with the exact nonlinear averages, 



. yjn\ + {2n h ) 2 a+a- ± n c yfi* + (2n b ) 2 e<M*)> + G [(S 3 (x))] ± n c 

(n±(x)) = = , (38) 

because of the nonvanishing quadratic and higher-order (y > 2) contributions of the electrostatic-potential 
deviations, 

5 v {x) = m-i>{x)] v . (39) 

In Appendix [A] we compare the uniform self-consistent linearized expansion densities (|3^) with the exact 
nonlinear averages (p8[). Justification of the neglect of the quadratic and higher-order contributions under 
linearized theory, which was done for the spherical case — but is trivially generalized for the planar case — 
is found in Appendix F of Ref. Q. However, we should mention that — although internal (within the semi- 
grand-canonical ensemble) self-consistency under linearization is achieved by using the uniform expansion 
densities ([56|) — global self-consistency (preserving the exact nature of the Legendre transformation between 
the semi-grand-canonical and canonical ensembles) requires also the inclusion of the quadratic contribution 
of the averages (^) , as discussed in detail in Appendix G of Ref. |jj . The linearized equations including 
these quadratic contributions to the expansion densities are presented in Appendix |d|, where it is shown that 
their inclusion do not improve the agreement between the linearized and the full nonlinear equations. 

Inserting the linearized equilibrium density profiles ( p7|) into the Poisson equation (||) yields the Debye- 
Hiickel-like (DH) equation,E9 

c\ 2 iI>(t) 9 

= n 2 mx) - ~V] + J [8(x + L) + S(x - L)] , (40) 

where the parameter 

n = — = ; =, (41) 

n++n_ ^nl + (2n b ) 2 

measures the relative importance of the counterions to the ionic strength in the interplane region \x\ < L, 

I=Un+ + n_) = \^nl + {2n h y = g = -7==. (42) 
2 2 2?7 y/i - jf 

The (effective) Debye screening length in the interplane region k _1 satisfies 

K 2 

k 2 = 8ir£ B I = b > k 2 , (43) 
V 1 -v 2 

showing that screening is enhanced compared to the salt reservoir. 

The gauge-invariant linearized electrostatic potential satisfying the DH-like equation ( ^0| ) subject to the 
boundary conditions (|l^) can be readily obtained, 

iK*) = +fi(l-*L , (44) 

\ sinn kL J 

where the average electrostatic potential for an arbitrary electrostatic surface potential ip^ is given by 

(V>(a:)) = ipL+rjKL&iKL), (45) 
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in terms of the Langevin function, 



Jzf (x) = cothx — — . 



(46) 



The linearized semi-grand-canonical potential, J7dh = ^DH["±(^)]oquii, is obtaining by inserting the 
equilibrium density profiles ( |37| ) — recalling that the self-consistent linearized averages, (n±(x)) — n±, 
are given by Eq. (36) — and the DH-like solution flli]) into the linearized semi-grand-canonical functional 
J? DH [n±(a:)], Eq. (35). After performing the integrations, we may cast the linearized semi-grand-canonical 
potential (per unit area) in the form 



DH 
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which yields the dimensionless excess linearized semi-grand-canonical potential per unit area, 
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written in terms of the dimensionless lengths 

k 1 = Kb K_1 j I = KbL, A = KbA. 

With these definitions we obtain the linearized self-energy o->dh(A, I — > oo) = 2/A 2 and 
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k 2 = k 2 (\,l) = v/1 + [2/(A/)] 2 = 



1 



^I + WIW y/l-V 2 (X,l) 

The dimensionless linearized osmotic-pressure difference is then given by 
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where we have made use of the total derivative, 



d d d?7 d dk d d rj , 2 \ ® ^V 2 9 
M = m + lil&v + ^M = di~~l 9r)~ ~2Tdk' 



and the derivative of the Langevin function, 



££'(x) = \- -^(x)- % 2 (x). 



(51) 



(52) 



(53) 



Alternatively, the same expression for the linearized osmotic-pressure difference may be obtained by 
performing a quadratic expansion of the gauge-invariant form of the nonlinear PB osmotic pressure ( |34| ) — 
similarly as obtained for the spherical case in Appendix E of Ref. 0, 



n DH = k 1 



l + ^i(Ci) + ^ 2 (0)-^ 2 (S 2 (x)) 



- 1, 



where the v-th order electrostatic-potential differences ( |39| ) read 

v ( cosh kx 
d v (x)=r] kL— — --1 
\ sinn k_L 



(54) 



(55) 



In the next section we will investigate the properties of the linearized osmotic-pressure difference defined by 
Eqs. (|U) or (fj) and compare it with its exact nonlinear counterpart (|33|). 
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4 Comparison of the exact nonlinear and the linearized equations 
of state 

As already pointed out in the literature, the linearized osmotic-pressure difference ITdh defined by 
Eqs. ( |5l"| ) or (|54|) yields artifacts in the low-temperature, large-separation or high-surface charge limits. 
In contradiction to the exact nonlinear result (p3|), which predicts that the osmotic-pressure difference is 
always positive, II > 0, the linearized version IlrjH becomes negative in the above mentioned limits. In an 
attempt to define the osmotic pressure in a linearized framework, Deserno and von GriinbergH introduced 
an additional (alternative) definition, III, cf. Eq. (43) of Ref. ||, that does not have the shortcoming of 
displaying any instabilities. On the other hand, we will show later that their partially unstable osmotic- 
pressure definition, cf. Eq. (44) of Ref. , coincides with the linearized version ( |5l| ) obtained in the previous 
section, II2 = Edh- Their general formulas, for the planar case (d = 1), need to be taken in the formal limit 
of vanishing volume fraction, <f> = a/ 1 — > — with a > being some arbitrary (reduced) length — which 
yields 
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1 = cosh'i/'opt 



sinh 2 y> opt 
2cosh^ op t 




Hi - 



2 

; 2 \smhkl 
PPi _ sinh 4 t/j opt ( 1 - a/l 
2n h 2 cosh 3 ^ opt «™ { 2a/l 
( k 2 l 2 



>0, 



Ifc2 4 

4 Vsinh 2 ^ 



kl coth kl 



V = K 1/2 (ka)I 1/2 (kl) - K 1/2 (kl)I 1/2 (ka), 
£ = K_ 1/2 (ka)I 1/2 (kl) + K 1/2 (kl)I_ 1/2 (ka), 



1 , £ 
— - + ka — 

V 2 V 



k 2 a 2 [ 1 



V 2 



(56) 



(57) 

(58) 
(59) 



where {/„, K u } are the modified Bessel functional^ of the first and the second kind, respectively, and i/) pt is 
the optimal linearization point, satisfying the relations 



tanh-0 op t = -r] 



coshV'opt = 



= k 2 



sinhi/'opt = 



2n h 



2 

A/' 



(60) 



In accordance with Eqs. (23) and (26) of Ref. |3), the two osmotic-pressure definitions can be recast in a 
simpler formal form, 



n x = k 2 
n 2 = k 2 



1 + 77^(0) + -6 2 (0) 



-1, 
1 



1 + 77^1(0) + ^ 2 (0)-^ 2 (S 2 (x)) 



(61) 
(62) 



from which one can see that the second osmotic-pressure definition coincides with the linearized osmotic- 
pressure difference ( j5l|) obtained in the last section, II2 = ilDH n while the first one III differs from Eq. (|Z|) by 
an omitted quadratic term. Analogously to the spherical casejil the term that distinguishes the two distinct 
osmotic-pressure definitions originates from the volume dependence of the optimal linearization point ipopt, 
as pointed out by Deserno and von Griinberg.0 

From its asymptotic expansions to be given next and its formal expression (|6l| ) , we see that III , although 
fully thermodynamically stable — related to its positiveness, cf. Eq. ([56|) — is inconsistent with a quadratic 
expansion of the gauge-invariant nonlinear PB pressure ([34]), because of the omitted last quadratic term 
of (|54|). Furthermore, we will show next that the consistent — although partially unstable — linearized 
osmotic-pressure difference II2 presents indeed a better agreement with the nonlinear osmotic pressure II in 
the weak-coupling and counterionic ideal-gas limits, when the underlying assumptions of the linearization 
are fulfilled. Therefore, although the alternative 111 displays the fortuitous advantage of preserving the 
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positiveness of the exact nonlinear pressure II, its derivation has no justification in our approach based on the 
minimization of the linearized semi-grand-canonical functional ^?dh [n± (x)] . Moreover, the partially unstable 
II2 corresponds indeed to the negative total derivative of the linearized semi-grand-canonical potential wdh 
with respect to the planes separation /, which we thus believe to be the consistent and correct definition of 
the osmotic pressure. 

Let us,now perform an explicit comparison between asymptotic expressions of II, the nonlinear osmotic 
pressureEl — obtained in Appendix [b| — and of the two corresponding linearized versions, ITi and II2, for 
the distinct regimes listed below. 

• Weak-coupling or zero-th order Donnan (Ib 0) limit: / — > 0, A — * 00, but finite product XI 
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(63) 
(64) 
(65) 



Counterionic ideal-gas limit: I — ► and finite A 



n = 


1 \ 

XI 


!- 


I 

3A + 


V45 


A 4 


iii = 


1 \ 

XI 


!- 


I 

3A + 


{-- 


A 4 


n 2 = 


-\ 

XI 


[- 


I 

5a + 


(-- 

V45 


A 4 



16 (l_ 

945 [x 

8 n 

189 VA 

8 n 

315 V A 



(//A) z 



(i/xy- 



(i/xy 



-1, 



(66) 
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Gouy-Chapman or high-surface charge limit: I — > and A/Z — > 
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(70) 
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Large-separation limit: I — > 00 and finite A 
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Looking at Eqs. © and @ one may see why the linearized osmotic-pressure difference n 2 becomes 
negative at the Gouy-Chapman and large-separation limits. In the Gouy-Chapman limit the leading term 
is given by the G (A~ 3 / 2 ), which is negative and overcomes the exponentially decaying G (A~ 2 ) term. The 
leading term of the large-separation limit is given by the G (l~ 3 ) contribution, which is negative and over- 
comes the three exponentially decaying lowest-order terms. In the full nonlinear solution, however, all 
algebraically decaying terms cancel in a nontrivial way, and eventually only an exponentially (positive) de- 
caying behavior is predicted. Note that both linearized versions, LT^ and II 2 , show asymptotic behaviours 
that disagree strongly from the nonlinear osmotic-pressure difference II. This clearly indicates that both 
linearized osmotic-pressure definitions are meaningless in these limits and so is the positiveness of III. 

We see that in the weak-coupling limit the self-consistent linearized osmotic pressure II2 and its nonlinear 
counterpart II agree up to the G (/ 4 ) terms, confirming the validity of the linearization when its underlying 
assumptions are fulfilled. The same occurs for the counterionic ideal-gas limit up to the G (I) terms. In both 
cases the fully stable IT has a worse agreement, one order lower than the partially unstable II2. However, in 
the large-separation limit, the two linearized and the nonlinear expressions disagree even qualitatively: the 
linearized asymptotics is algebraic (negative for II2, positive for IT), whereas the nonlinear is exponential 
(and positive). On the other hand, although in the Gouy-Chapman limit all asymptotics are algebraic, 
in the linearized case the power-law is oc — l~ x / 2 for II2 and oc l~ x for III, both in disagreement with the 
nonlinear asymptotics cx l~ 2 . The failure of the linearization scheme should not be at all surprising, because 
it is supposed to be valid in the weak-coupling (£b — * 0) and counterionic ideal-gas limit, but not in the 
opposite, large-separation (I — > 00) or high-surface charge, Gouy-Chapman, (A — > 0) limits. Therefore, any 
results obtained in a linearized framework outside the weak-coupling and the counterionic ideal-gas limits 
should be taken with caution. See also our comments in the concluding remarks of the preceding papeiu on 
the gas-liquid-like phase separation in dilute deionized aqueous suspensions of charged colloidal particles. 

In order to show the accuracy of the self-consistent linearized osmotic-pressure difference IIdh , Eq. ( |5l| ) , 
and the region where the linearization scheme breaks down, we plotted in Figures [l] and || the locii of constant 
errors between the exact nonlinear PB osmotic-pressure difference and the corresponding linearized version, 
measured by the logarithmic deviations 

OTdh = |lnn DH (A,Z)-lnn(A 1 0|. (75) 

We have chosen a logarithmic measure for the deviations because II varies in a range of several orders of 
magnitude. For small deviations, this definition leads to the relative errors, 



sn 



DH 



n DH (A,0-n(A,0 



n(A,o 



(76) 



Analogously, we may define the logarithmic deviation from PB of the linearized semi-grand-canonical po- 
tential, 

dujun = |ln[u/DH(A,0 - w(A,oo)] - ln[w(A,/) - w(A,oo)]| , (77) 

which is always smaller than 5IIdh (not shown). Therefore the linearized semi-grand-canonical potential, 
Eq. ([f8|), and the linearized osmotic-pressure difference, Eq. ([5l]) , describe well the corresponding nonlinear 
equations in the limit A/I > 1 and I « 1. Because the nonlinear theory always predicts repulsion, the 
attractive osmotic-pressure region — shown in gray in Figures Q and || — is clearly an artifact of the 
linearization. When plotted on the A/L x (k^L)" 1 plane, the IIdh = line reaches at k^L — ► the 
asymptotic value £0 = A/L — 0.123863965 • • •, which is given by the solution of the transcendental equation 

%o + Mj?(\[f) +^ 2 (\f^) =1- (78) 



2 VVW VV Co 

To obtain the full nonlinear PB osmotic-pressure difference II, one needs to numerically solve the tran- 
scendental equation (|2^) involving elliptic functions or elliptic integrals. Although the asymptotic expansions 
of the nonlinear II represented by Eqs. (|63|), (|66|), (|69| ) and (|72|) allow an explicit analytical comparison in 
the distinct regimes with their linearized versions, they are not very useful for numerical evaluation. In Ap- 
pendix |c] we derive extended expansions of the nonlinear PB semi-grand-canonical potential u>, Eq. (|30|), and 
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of the PB osmotic-pressure difference II, Eq. (p3|), that involve only elementary functions and are suitable 
for numerical implementation. These extend the numerical accuracy of the above mentioned asymptotic 
expansions of the full nonlinear II and are complementary to the linearized equations, udh, Eq. (fl8|), and 
IIdh, Eq. (|5l|), providing an excellent approximation in the regions where the linearization scheme breaks 
down. In Figures |l| and p| we also present their corresponding logarithmic deviations from the exact PB 
result, which, similarly to (|75[), are defined by 

5U GC = |mn GC (A,0-mII(A,0|, (79) 
Sn hs = |lnn LS (A,0-lnII(A,0|, (80) 

where Hgc and Bils, given explicitly in Appendix |^, are the osmotic-pressure differences in the extended 
Gouy-Chapman and large-separation limits, respectively. 



5 Concluding remarks 

The classical problem of two infinite uniformly charged planes in electrochemical equilibrium with an infinite 
salt-reservoir is exactly solved at the mean-field nonlinear level, as well as by a linearization scheme consistent 
with quadratic expansions of the nonlinear semi-grand-canonical functional. By using gauge-invariant forms 
of the electrostatic potential, we have shown that the linearized osmotic pressure corresponds to a quadratic 
expansion of the corresponding nonlinear version. 

As already pointed out in the literaturejj it is shown that the self-consistent linearized osmotic pressure 
leads to artifacts in the large-separation and the Gouy-Chapman (high-surface charge) limits, predicting 
there negative osmotic-pressure differences. Although it is possible to define an alternative linearized osmotic 
pressure that it is fully stable based on the partial derivative of the linearized semi-grand-canonical potential 
with respect to the separation distancep its stability is shown to be a fortuitous result. In fact explicit 
comparison of the exact nonlinear osmotic pressure and the two linearized versions shows that the linearized 
self-consistent osmotic pressure, though partially unstable, presents a better agreement with the PB results 
in the weak-coupling and counterion ideal-gas limits, where the linearization can be applied. However, not 
surprisingly, in the region where the linearization breaks down none of both proposed linearized osmotic 
pressures give quantitatively correct results. 

To avoid confusion we should stress at this point the exactness of the PB nonlinear solution at the mean- 
field level and discuss its range of validity and limitations. It is known from numerical simulations of the Prim- 
itive ModeE3 (PM) in the planar geometry that sufficiently close and highly charged planes in the presence of 
neutralizing counterions attract each other ,L3 even though for realistic charge densities and monovalent ions 
this is not observed at room temperature. In this case the attraction is prevented by steric repulsions at the 
small separations at which it would be observed neglecting the finite ionic size. Because the mean-field PB 
approximation always predicts repulsion, theoretical validation for this attraction (observed in fact at room 
temperature only for multivalent ions) has_to be given beyond the PB level, e.g., b^hulk counterion corre- 
lations,tJ'H integralp£imatipns theories ,Oc3 charge-correlation-inducecLattractions,L3li3 charge-fluctuation- 
induced attractions Jlj'OcZI electrolytic depletion-induced attractions ,C3 discrete solvent-mediated attrac- 
tions,^! field-theory methodsElJ etc — see also Refs. H,H,||||,||,||,||,||||,|| for mechanisms of 
attraction between like-charged rods. On the other hand, in the strong-coupling limit the linearization of 
the WS-cell mean-field PB equation, as discussed in this work, does predict attraction. However, here the 
mechanism of attraction is related to mathematical artifacts of the linearization itself and does not corre- 
spond to a real physical effect. The fact that this prediction is in agreement with the theories beyond the 
mean-field level is purely accidental and is intrinsically connected with the inadequacy (meaning incorrect 
application) of the PB mean-field approach at the same limit. In other words, a qualitatively correct result 
(in this example, attraction) may be deceptively anticipated in the strong-coupling limit because of the 
simultaneous application of two inadequate approximations, namely, the mean-field PB equation and its 
subsequent linearization. 
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A Exact nonlinear averaged densities 



In this Appendix we will compare the uniform expansion densities about which the linearization is performed 
— the state-independent zero-th order Donnan densities — with the exact nonlinear PB averages. 
By using the definite integrals 
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it is possible to obtain the exact nonlinear PB averaged densities, 
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In Figure |jj we compare them with the uniform densities about which the linearization is performed, the 
state-independent zero-th order Donnan densities (pf 



(A4) 
(A5) 



by looking at their logarithmic deviations from the corresponding exact PB averages, 

5n± = In (n± (x)) — In n± . 



B Asymptotic expansions of the nonlinear solution 



In this Appendix we obtain the asymptotic expansions of the nonlinear osmotic-pressure difference II. We 
have made extensive use of Refs. 18, 19[^o[^l[ throughout this Appendix. 



B.l Weak-coupling limit 

To obtain the weak-coupling (Zb — ► 0) limit, first note that the product XI — 2/(rjk 2 ) does not depend on 
Ib- Therefore we multiply both sides of the eigenvalue equation by I, 



XlVt = ^(l = I cs ( ZVT+I 
r\k z 



1 



l + t 



(Bl) 



and expand them in powers of I, assuming that t = T^kLo a 2kl 2k - It is then possible to obtain the coefficients 
{o-2k] of the expansion recursively, leading to Eq. 



B.2 Counterionic ideal-gas limit 

To obtain the counterion-dominated ideal-gas limit (I — > and finite A), we write 

1 



XVt = cs Wl + t 



l + t 



(B2) 



and expand the right-hand side in powers of I, assuming that t = y^^° = _ 1 a,kl k ■ It is then possible to obtain 
the coefficients {a^} of the expansion recursively, leading to Eq. (p6[). 
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B.3 Gouy-Chapman limit 

To obtain the Gouy-Chapman (high-surface charge) limit (A — > 0), note that 
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If we additionally assume that / — > 0, with X/l — > 0, and with the help of the expansion of the complete 
elliptic integral K{m) about m = 0, 
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l + ^ + ^(m 2 )' 



we are lead to the asymptotic solution 
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which gives the asymptotic osmotic-pressure difference in the Gouy-Chapman limit, Eq. (p9|). Evaluation of 
higher-order terms of the leading contribution ff(l~ 2 ) of the expansion would give 
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B.4 Large-separation limit 

To obtain the large-separation (I — > oo and finite A) asymptotics, first note that 
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Using the expansion of the complete elliptic integral K(m) about m = 1, 
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we are led to the asymptotic solution 
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which yields the leading term of the asymptotic osmotic- pressure difference in the large-separation limit, 
Eq. (|72|). By computing higher-order terms, it is possible to write an asymptotic power series, 
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where the coefficients t u (X, I) = (D(l u ) of the three leading terms read 
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Although this asymptotic series works pretty well for large separations (I ^S> 1), in the crossover region (I w 1) 
it leads to oscillating pressures. In Appendix [c] we obtain extended expansions of the nonlinear equations 
in the large-separation limit that do not have this disadvantage in the crossover region. 



C Extended expansions of the nonlinear solution 

In this Appendix we obtain extended expansions of the nonlinear semi-grand-canonical potential and of the 
nonlinear osmotic-pressure difference that are valid in the region where the linearization breaks down. Again, 
We have made extensive use of Refs. [18ll9^,EG,El| throughout this Appendix. 



C.l Extended Gouy-Chapman limit 

In the previous Appendix, both the counterionic ideal-gas (finite A) as well as the Gouy-Chapman (A — ► 0) 
asymptotics were obtained in the small-separation (I — > 0) limit. In fact, for any ratio £ = X/l the summation 
over the X/l series for the leading terms up to 0{l 2 ) may be performed exactly, yielding 
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where y — y(£) is the solution of the transcendental equation 

&/tany = 1. 



l 2 + ff^), 



(CI) 



(C2) 



This general expression yields the leading term 0(l~ 2 ) of the counterionic ideal-gas (finite A, when y — > 
\fl/X — > 0) as well as the Gouy-Chapman (high-surface charge, when A — > 0, y — > tt/2) asymptotics as 
special cases, 
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The excess semi-grand-canonical potential ui may be obtained by integration of the osmotic-pressure differ- 
ence 2t, leading to 
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While the third and fourth terms are the leading corrections due to the presence of salt, the two first terms 
can be related to half of the exact nonlinear Helmholtz free energy of two charged infinite planes in the 
presence of neutralizing counterions only (salt-free Gouy-Chapman case), 
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where y is the solution of the transcendental equation y tany = L/A. 

We define the extended Gouy-Chapman limit by truncating the above expansions, neglecting thus higher- 
order terms, 
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where y = y(£) is the solution of the transcendental equation (|C2|). 
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C.2 Extended large-separation limit 



The large-separation osmotic-pressure asymptotics (Bll) obtained in Appendix B.4 displays oscillations in 
the crossover (I w 1) region. Because we want to match the linearized DH-like, the extended Gouy-Chapman 
and the large-separation asymptotic expressions at the crossover region, we need to find an extended expan- 
sion that does not display this shortcoming. In fact the pressure oscillations are avoided if one uses instead 
the implicit form I = /(A,m), which is obtained by expanding the eigenvalue equation, 
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in powers of (1 — to) = t/(l + t). Accurate results in the crossover region, which will cover almost the whole 
(/ x A) parameter space with logarithmic pressure deviatipns from the exact PB within 0.1%, are obtained 
by using fourth-order expansions of the elliptic integrals^ about to = 1, 
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The asymptotic large-separation (Z — > oo) excess semi-grand-canonical potential is obtained by integration 
of the osmotic-pressure difference, 2t = — div/dl, 
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(C12) 



(C13) 



which can all be performed analytically. After some algebraic manipulations, the last contribution of 
Eq. (C12) may be cast in the form 



l M 



dfj, 315 
— l(X,fx) = 



32 



arctanh y/m + In 2 — In ( A 



ln(l-m) -41n2 + 21n(A+ Vl + A 2 
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56 



C(m) 



XA(X) 



B{m) 



3360(1 + A 2 ) 7/2 8192^ 
^A 2 " +1 C„(to) 



49152^ 430080 (1 + A 



2^/2 



(C14) 



A(X 
B(m 
C(m 
Co (to 
Ci(to 
C 2 (m 
C 3 (m 
Ca(jti 
C 5 (to 
C 6 (to 
C?(m 



34125 + 24150A 2 - 56210A 4 - 46480A 6 + 9253A 8 + 738A 10 - 5160A 12 + 1680A 14 , (C15) 



= 25609 + 18404to - 4818to 2 + 1300to 3 - 175to 4 , 

= 123743 + 82792to - 56970to 2 + 16960to 3 - 2365to 4 , 

= 105 (14711 + 35356to - 11310to 2 + 3308to 3 - 465m 4 ) , 

= 70 (19795 + 26252to - 1254to 2 - 868m 3 + 235m 4 ) , 

= -70 (31873 + 98756m - 37826m 2 + 11636m 3 - 1655m 4 ) , 

= -560 (4129 + 7940m - 1730m 2 + 308m 3 - 23m 4 ) , 

= 156695 + 1690780m - 938350m 2 + 322924m 3 - 47665m 4 , 

= 79030 - 82600?7i + 186900m 2 - 111496m 3 + 22630m 4 , 

= -40 (4795 + 16940m - 7350m 2 + 2492m 3 - 365m 4 ) , 

= 1680 (35 + 140m - 70m 2 + 28m 3 - 5m 4 ) . 



(C16) 
(C17) 
(C18) 
(C19) 
(C20) 
(C21) 
(C22) 
(C23) 
(C24) 
(C25) 



We may check numerically that the osmotic pressure 2t corresponds indeed to the derivative of the excess 
grand potential lu with respect to the semi-separation I between the two charged surfaces. 

The extended large-separation expressions are obtained by truncating the above expressions, neglecting 
thus higher-order contributions. Collecting all contributions, the extended large-separation osmotic-pressure 
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difference IIls = 2(1 — to) /to is defined implicitly by the relation 



l(X, to) 



ln4 ln(l — to) H — (1 — to) 

2 4 



ln4 — 1 — — ln(l — m) 



+ — (1 - m) 2 [6 In 4 - 7 - 31n(l - to)] + t^t(1 - m) 3 [30 In 4- 37- 151n(l - m)] 
128 1536 



35 



196608 
1 



(1 - m) 4 [420 In 4 - 533 - 2101n(l - to)] - 



A (1 - to) 



2(1 + A 2 ) 



h{( i+x2 r 



(1 - to) (1 + A 2 ) 2 (3 - 2A 4 ) + (1 - to) 2 (1 + A 2 ) (15 + 5A 2 - 10A 4 + 6A 6 + 9A 8 ) 



1 
192 



(1 - to) 3 (105 + 70A 2 - 70A 4 + 28A 6 + 66A 8 - 72A 



10 



60A 



12> 



In 



(A + Vl + A 2 



1 — TO 



ln(A+ Vl + A 2 ) -AVl + A 2 



31n(A+ VT+A 2 ) - AVl + A 2 (3 - 2A 2 



15 In (A + Vl + A 2 J - AVTTA 2 (15 - 10A 2 + 8A 4 
105 In (A + Vl + A 2 ) - A^l + A 2 (105 - 70A 2 + 56A 4 - 48A 6 



(C26) 



(C27) 



3(1 -to) 
64 

5 (1 - mf 

768 
35(1 - to)' 

49152 

with the associated semi-grand-canonical potential, 

WLs(A,0s|arccosh^l + ^ + 4 (\ - j^l + A 2 ^ - |- (1 - to) Z(A, to) - 2 ^ ^l(X,fx). 

D Globally self-consistent linearized equations 

In this Appendix we show that the linearized equations that preserve the exact nature of the Legendre 
transformation do not lead to any improvements in the agreement between the linearized and nonlinear 
osmotic pressures in comparison to the linearized versions obtained in Section ||. 

As discussed in detail in Appendix G of Ref. [Q , the Legendre transformation between the semi-grand- 
canonical and the canonical descriptions of the system may be rendered exact if — instead of using the 
state-independent zero-th order Donnan densities ( |36| ) — one uses the quadratic truncation of the nonlinear 
averages, given by Eq. (|38|), as expansion densities to obtain the linearized functional. With the inclusion 
of the quadratic state-dependent contribution (^(x)) to the average densities, we obtain the globally self- 
consistent linearized semi-grand-canonical potential and linearized osmotic-pressure difference, 



cj D h(A,0 = 



A 



, ~ 1 5 
arctanh 77 1 — 

77 4 



2n^L 



1 1 



arctanh r) — + — fjklJ? I kl\ -\ (J 2 (x)) 



1 



2f) 



+ 1 



-fi(kl) 2 



j£f 2 (kl) - 1 



+ 1, 



n DH (A,/) 



dw DH (A,0 _ ~ 2 



d/ 



1 + ^r(Q) + I<5 2 (0) - i (5 2 (x)) 



pa 



^ffk& (w) + i(#o 2 ^ 2 (w) 



- 1 



- 1, 



where the dimensionless parameters, 

1 



v /l + (A//2) 2 e^( :E ))' 



fc 2 = y/ e (B a (x)) + [2/(AZ)] 2 = 



e <«2(x)>/2 



(Dl) 



(D2) 



(D3) 
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are given implicitly in terms of the quadratic electrostatic-potential deviation, 



(5 2 (x)) = ^fk (kij + k ^f 2 (w) - 1 } 



(D4) 



To compute the linearized osmotic-pressure difference IIdh, Eq. (D2), one needs to take into account the 
total derivatives of the parametric forms, Eqs. (D3), 



d d dfi d dk d 
_ = 1 1 

dl dl dl dr) dl dk 
d 



^ Ci - 



c) 



1 , I d(M) 

2 dl J df) 



~Zl 



1 _ I d <w> 



d_ 

dk 



(D5) 



In accordance to the infinite-separation linearized self-energy obtained in Section ^ the new version is also 
given by w D h(A, I ^ 00) = 2/A 2 . 

In Figures Fjjand || we compare the two linearized osmotic-pressure definitions, IIdh and IIdh, given by 
Eqs. ( [si| ) and (|P2|), with the exact nonlinear version II, given by Eq. (^3|). The dotted lines in Figures [l] and || 
suggest a better agreement between the linearized osmotic-pressure difference IIdh and the full nonlinear 
counterpart II — in comparison to the linearized version IIdh, obtained in Section |[ However, as shown 
below by explicit analytical comparison, these numerical evidences are in fact misleading. 

Asymptotic analytical expansions of the linearized osmotic-pressure difference IIdh about the weak- 
coupling (£ B -> 0), 



n 



DH 



IIdh + k 

,2 



16200 



(r? 2 + 5) (l- v 2 )k 8 l 8 + ff(k 10 l 10 ) 



6 90 v ' ; 



945 v ' ; 113400 



(7?7 6 + 28t7 4 + if 



= k 2 
k»l* + ff(k w l w ) 



v t 2 



and the counterionic ideal-gas (L — > 0, finite A) limits, 

IIdh = IIdh 




(i/xy 



(D6) 



(D7) 



show explicitly that both linearized osmotic pressures, IIdh and IIdh, agree with the full nonlinear PB 
version II up to the same order — cf. Eqs. (|63|) to (|68|). Therefore, the numerical indications of a better 
agreement of IIdh, as suggested by Figures |l] and |[ are purely fortuitous. In fact, for ratios A/L > 10 2 
(beyond the values shown in Figure ||) one observes a crossover between the deviations of the linearized 
versions, IIdh and IIdh, with respect to the full nonlinear osmotic- pressure difference II. 
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Figure 1: Logarithmic deviations from the PB of the different asymptotic osmotic-pressure differences. 
DH represents the Debye-Hiickel-like, linearized-functional expansion about the weak-coupling limit, GC 
corresponds to the expansion about the salt-free Gouy-Chapman limit, and LS denotes the large-separation 
limit expansion. The region complementary to LS is splitted into three parts. In the gray region the 
linearized osmotic-pressure difference IIdh becomes negative. The arrows indicate the direction of decreasing 
logarithmic deviation 5H from the PB results: 1CP 1 (dot-dashed lines), 10~ 2 (dashed lines), 10~ 3 (solid lines). 
For comparison, we also display (dotted lines) the linearized results by i nclu ding quadratic contributions in 
the expansion densities, as defined by the linearized pressure LIdHj Eq. ( |D2| ). 




1 2 3 4 5 6 

Figure 2: Same as in Figure [j], but plotted using different variables. At n^L — > 0, the IIdh = line reaches 
the asymptotic value defined by Eq. (fTSt), £o = A/L = 0.123863965 • • • Compare with Figure 1 from Ref. p5| . 
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l=K b L 

Figure 3: Deviations from the PB averaged densities of the state- independent zero-th order Donnan densities, 
which were used to perform the quadratic expansions of the nonlinear functional. The arrows indicate the 
direction of increasing logarithmic deviations 5n± from the PB results: 1CP 3 (solid lines), 10 -2 (dashed 
lines) and 10 _1 (dot-dashed lines). To allow a comparison with the region where the linearized theory 
breaks down, we also plotted the locus (dark gray thick line) where the logarithmic deviation from PB of 
the linearized osmotic-pressure difference is <5I1dh = 1CP 3 . In the light gray region the linearized osmotic- 
pressure difference becomes negative, PIdh < 0. Although there is a close connection between this region 
and the increase of the deviations 5h± for high-surface charges (A <C 1), for low-surface charges (A 3> 1) and 
large separations (I ^> 1), the linearized theory still predicts a negative linearized osmotic-pressure difference 
(upper-right region), while the full nonlinear one vanishes exponentially from positive values. 
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